#include "PowerSpectrum.h"
#include "Int.h"

#include <fstream>
using std::ifstream;
using std::ofstream;
#include <iostream>
#include <sstream>
using std::istringstream;
#include <string>
using std::string;
#include <map>
using std::map;
#include <vector>
using std::vector;

#include "dbg.h"
using std::endl;

//#define STARTAT 68955

std::ostream* dbgout=0;
//bool XDEBUG=true;
bool XDEBUG=false;

//const char* cmbfile = "CMBQuintchainNew.dat";
//const char* cmbfile = "wmap3_wcdm_pert.dat";
//const char* outfile = "pred.wmap3.out";
//const char* debugfile = "pred.wmap3.debug";

#define WMAP3

//const char* cmbfile = "WWaPrior.dat";
//const char* outfile = "pred.wwa.out";
//const char* debugfile = "pred.wwa.debug";

#define LAMBDA
//const char* cmbfile = "LCDMWMAP.dat";
const char* cmbfile = "wmap3_lcdm.dat";
const char* outfile = "pred.lambda.out";
const char* debugfile = "pred.lambda.debug";

//const char* cmbfile = "CMBQuintchain.dat";
//const char* outfile = "pred.w1.out";
//const char* debugfile = "pred.w1.debug";

//#define Wa
//#define NonFlat
//#define TAKEEVERYN 100 // Only take every nth line from cmbfile

//const char* cmbfile = "NonFlatWWaPrior.dat";
//const char* outfile = "pred.nfwwa.out";
//const char* debugfile = "pred.nfwwa.debug";

const double PI = 3.141592653589793;
const double DEG = PI/180.;
const double ARCSEC = DEG/3600.;
const double ARCMIN = DEG/60.;

const double map_minr = 1.0*ARCMIN;
const double map_maxr = 40.0*ARCMIN;
const double var_minr = 10.0*ARCMIN;
const double var_maxr = 100.0*ARCMIN;
const double stepfactor = 2.;

#ifdef NonFlat 
typedef map<UParamsNonFlat,double> umap;
#else
typedef map<UParams,double> umap;
#endif

void ReadCMB(umap& lnlike)
{
  dbg<<"Start ReadCMB\n";

  dbg<<"Using "<<cmbfile<<" for priors\n";
  ifstream fin(cmbfile);
  if (!fin) myerror("opening cmb file:  ",cmbfile);

  double ombh2,omch2,theta,n,Z,w=-1.,As,h,like;

  // I ignore theta and Z. 
  // ll is apparently not acually the log likelihood.
  // Licia says to ignore it and take the likelihood
  // to be proportional to the repetition of the lines.
 
#ifdef WMAP3
  double junk;
  //                   -ln(lik)                  theta    tau     om_k
  while (fin >> like >> junk >> ombh2 >> omch2 >> junk >> junk >> junk >>
  //  f_nu              n_t     dn_d/dlnk      r      A_sz    om_lam
      junk >> w >> n >> junk >> junk >> As >> junk >> junk >> junk >>
  //  age     om_m    sig_8   z_rei   ZERO       beta_sdss b_sdss  z_bg
      junk >> junk >> junk >> junk >> junk >> h >> junk >> junk >> junk)
#else
  while (fin >> ombh2 >> omch2 >> theta >> n >> Z >>
#ifndef LAMBDA
      w >> 
#endif
      As >> h >> like) 
#endif // WMAP3
  {

#ifdef WMAP3
    h /= 100;
    // The input As is really ln(1.e10 * 2.95e-9 * As)
    // = ln(29.5 * As)
    As = exp(As);
    As /= 29.5;
    //dbg<<"h,As = "<<h<<", "<<As<<endl;
#endif

#ifdef Wa
    double wa = Z;
#else
    double wa = 0.;
#ifndef WMAP3
    like = 1.;  
    // For Wa, I calculate likelhood for each line
    // And for WMAP3, like is the multiplicity of lines
#endif
#endif
    
#ifdef NonFlat
    double omde = theta;
    UParamsNonFlat up(ombh2,omch2,omde,n,w,wa,As,h);
#else
    UParams up(ombh2,omch2,n,w,wa,As,h);
#endif

    //xdbg<<"om_m = "<<up.GetOmM()<<", w = "<<w<<", wa = "<<wa<<endl;
    lnlike[up]+=like;  // for now, _not_ log
  }
  dbg<<"cmblnlike size = "<<lnlike.size()<<endl;
  umap::iterator mapit;
  for(mapit = lnlike.begin(); mapit != lnlike.end(); mapit++) {
    mapit->second = log(mapit->second);
    //dbg<<mapit->first.w<<"  "<<mapit->first.As<<"  "<<mapit->second<<endl;
  }
}

int main(int argc, char* argv[]) 
{
  dbgout = new ofstream(debugfile);
  //dbgout = &std::cout;

  dbg<<"map range = "<<map_minr/ARCMIN<<" - "<<map_maxr/ARCMIN<<endl;
  dbg<<"var range = "<<var_minr/ARCMIN<<" - "<<var_maxr/ARCMIN<<endl;

  umap cmblnlike;
  ReadCMB(cmblnlike);

  umap::iterator mapit;
  ofstream fout(outfile);
  int nmap = int(floor(log(map_maxr/map_minr)/log(stepfactor))) + 1;
  dbg<<"nmap = "<<nmap<<endl;
  double mapstep = exp(log(map_maxr/map_minr)/nmap);
  dbg<<"mapstep = "<<mapstep<<endl;
  int nvar = int(floor(log(var_maxr/var_minr)/log(stepfactor))) + 1;
  dbg<<"nvar = "<<nvar<<endl;
  double varstep = exp(log(var_maxr/var_minr)/nvar);
  dbg<<"varstep = "<<varstep<<endl;

  fout << map_minr<<"  "<<map_maxr<<"  "<<nmap<<"  "<<mapstep<<endl;
  fout << var_minr<<"  "<<var_maxr<<"  "<<nvar<<"  "<<varstep<<endl;

  int count=0;
  for(mapit=cmblnlike.begin();mapit!=cmblnlike.end();mapit++,count++) {

#ifdef STARTAT
    if (count < STARTAT) continue;
#endif
#ifdef TAKEEVERYN
    if (count % TAKEEVERYN != 0) continue;
#endif

#ifdef NonFlat
    const UParamsNonFlat& up = mapit->first;
#else
    const UParams& up = mapit->first;
#endif

    //double cmbll = mapit->second;
    dbg<<count<<"/"<<cmblnlike.size()<<"  "<<up.ombh2<<"  "<<up.omch2;
#ifdef NonFlat
    dbg<<"  "<<up.omde;
#endif
    dbg<<"  "<<up.n<<"  "<<up.w<<"  "<<up.wa<<"  "<<up.As<<"  "<<up.h<<endl;

    Universe u(up);

    double sig8 = u.sigma8();
    dbg <<"sigma8 = "<<sig8<<"  Omega_m = "<<u.om_m<<endl;

    bool settozero = false;
    double sig8_03 = sig8 * sqrt(u.om_m/0.3);
    dbg<<"At Omega_m = 0.3: sig8 = "<<sig8_03<<endl;
    if (sig8_03 < 0.2 || sig8_03 > 2.6) {
      dbg<<"Set map and var predictions to all 0's\n";
      settozero = true;
    } else if (sig8_03 < 0.4 || sig8_03 > 1.3) {
      dbg<<"Sloppy calc is ok.\n";
      u.precise = false;
    }

    double r = map_minr;
    for(int i=0;i<=nmap;i++, r*=mapstep) {
      double map = settozero ? 0. : u.fmm(r);
      fout << map <<"  ";
      dbg<<"map("<<r/ARCMIN<<") = "<<map<<endl;
    }
    r = var_minr;
    for(int i=0;i<=nvar;i++, r*=varstep) {
      double var = settozero ? 0. : u.fvar(r);
      fout << var <<"  ";
      dbg<<"var("<<r/ARCMIN<<") = "<<var<<endl;
    }
    fout<<sig8<<endl;
  }

  return 0;
}
